Este es el pipeline seguido para el análisis de los metagenomas de amplicones de ITS del trabajo de Daniella Senatore, utilizando dada2. Para correrlo se requieren los paquetes dada2, ShortRead y phyloseq, así como también Qiime2. Para instalar Qiime2 podemos consultar el tutorial disponible en https://docs.qiime2.org/2020.11/install/
dada2 requiere que los archivos que contienen las secuencias sean previamente desmultiplexadas. Para esto, hay diferentes opciones disponibles. Yo utilicé Qiime2, siguiendo los siguientes pasos (se asume que Qiime2 está instalado y accesible):
Antes de importar las secuencias a Qiime debemos comprimir el archivo utilizando el programa gzip:
Requiere:
Se obtiene: archivo .fastq.gz
gzip archivo_fastq_multiplex.fastq
Dado que utilizaremos Qiime2 para demultiplexar los reads, el primer paso en el pipeline de Qiime es la importación.
Requiere:
Se obtiene:
qiime tools import --type MultiplexedSingleEndBarcodeInSequence --input-path archivo_fastq_multiplex.fastq.gz --output-path multiplexed-seqs.qza
Argumentos del comando:
--input-path: path al archivo que se importa.
--output-path: nombre que se le da al archivo importado. Tiene extensión qza.
--type MultiplexedSingleEndBarcodeInSequence: indica que el archivo que se importa es multiplexado, contiene barcodes y contiene reads single end (SE).
Output:
Saved SampleData[SequencesWithQuality] to: demultiplex_seqs.qza
Saved MultiplexedSingleEndBarcodeInSequence to: unmatched_barcode.qza
Este paso se realizó según las instrucciones disponibles en https://forum.qiime2.org/t/demultiplexing-and-trimming-adapters-from-reads-with-q2-cutadapt/2313.
Requiere:
Se obtiene:
qiime cutadapt demux-single --i-seqs multiplexed-seqs.qza --m-barcodes-file mappingfile.txt --m-barcodes-column BarcodeSequence --o-per-sample-sequences demultiplex_seqs.qza --o-untrimmed-sequences unmatched_barcode.qza --output-dir demultiplex
Argumentos del comando:
cutadapt demux-single: se ejecuta el programa Cutadapt para demultiplexar secuencias SE.
--o-untrimmed-sequences: locación de las secuencias que no matchearon con ningún barcode.
--i-seqs: nombre del archivo importado (qza), que se generó en el paso anterior.
--m-barcodes-file: path al mapping file.
--m-barcodes-column: nombre de la columna del mapping file que contiene la secuencias de los barcodes.
Una vez demultiplexados, deben removerse los barcodes de los reads. Para esto volveremos a utilizar el plug-in Cutadapt de Qiime.
Requiere:
Se obtiene:
qiime cutadapt trim-single --i-demultiplexed-sequences demultiplex_seqs.qza --o-trimmed-sequences trimmed-demult_seqs.qza
Argumentos del comando:
cutadapt trim-single: se ejecuta el programa Cutadapt para trimmear secuencias de adaptadores o barcodes en reads SE.
--i-demultiplexed-sequences: archivo .qza que contiene los reads demultiplexados obtenidos en el paso anterior.
--o-trimmed-sequences: nombre del archivo de salida (qza) que contiene los reads demultiplexados y trimmeados.
Hasta ahora venimos trabajando con Qiime2, por lo tanto todos los archivos que se generan durante el procesamiento de las muestras está almacenado en artefactos. Para continuar el procesamiento de los datos en dada2 debemos transformarlos en fastq. Para esto exportamos los datos contenidos en archivos .qza.
Requiere:
Se obtiene:
mkdir demulti
qiime tools extract --input-path trimmed-demult_seqs.qza --output-path trim-demult
Argumentos:
--input-path: path del archivo qza a exportar
--output-path: nombre de la carpeta que contiene los datos exportados
dada2En este punto contamos con los reads desmultiplexados y sin barcodes, y continuaremos su procesamiento en dada2. Si los paquetes necesarios no están instalados, pueden obtenerse ejecutando:
# Librerías auxiliares
library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
── Attaching packages ────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
✓ ggplot2 3.3.2 ✓ purrr 0.3.4
✓ tibble 3.0.4 ✓ dplyr 1.0.2
✓ tidyr 1.1.2 ✓ stringr 1.4.0
✓ readr 1.4.0 ✓ forcats 0.5.0
── Conflicts ───────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::collapse() masks Biostrings::collapse(), IRanges::collapse()
x dplyr::combine() masks Biobase::combine(), BiocGenerics::combine()
x purrr::compact() masks XVector::compact()
x purrr::compose() masks ShortRead::compose()
x dplyr::count() masks matrixStats::count()
x dplyr::desc() masks IRanges::desc()
x tidyr::expand() masks S4Vectors::expand()
x dplyr::filter() masks stats::filter()
x dplyr::first() masks GenomicAlignments::first(), S4Vectors::first()
x dplyr::id() masks ShortRead::id()
x dplyr::lag() masks stats::lag()
x dplyr::last() masks GenomicAlignments::last()
x ggplot2::Position() masks BiocGenerics::Position(), base::Position()
x purrr::reduce() masks GenomicRanges::reduce(), IRanges::reduce()
x dplyr::rename() masks S4Vectors::rename()
x dplyr::slice() masks XVector::slice(), IRanges::slice()
x tibble::view() masks ShortRead::view()
Una vez instalados los paquetes requeridos se cargan en la sesión de R, utilizando la función library():
library(dada2)
Loading required package: Rcpp
library(ShortRead)
library(phyloseq)
En el paso anterior exportamos los datos procesados por Qiime. Ahora guardamos en una variable llamada path la ubicación de los archivos exportados. En este caso, los archivos se encuenctran en el directorio trimm_demult:
path <- "trimm_demult"
Ahora creamos el objeto fnFs, que contiene los path a los archivos fastq demultiplexados y sin barcodes, especificando un patrón que tengan en común los nombres de los archivos. En este caso ese patrón es __001.fastq.gz, ya que todos los archivos a analizar contienen este sufijo.
fnFs <- sort(list.files(path, pattern = "_001.fastq.gz", full.names = TRUE))
Los reads procesados anteriormente aún contienen las secuencias adaptadoras que se adicionaron a los fragmentos de ADN durante la creación de las librerías para la secuenciación, que debemos conocer de antemano. Además, para que dada2 funcione no puede haber ningúna base indeterminada (N) en las secuencias.
Primero, creamos un directorio donde vamos a guardar los reads sin Ns, y guardamos la ruta a los archivos que se crearán en un objeto:
dir.create("filtN")
fnFs.filtN <- file.path(path, "filtN", basename(fnFs))
fnFs.filtN
[1] "trimm_demult/filtN/A_AGGCAATTGC_L001_R1_001.fastq.gz"
[2] "trimm_demult/filtN/B_TTAGTCGGAC_L001_R1_001.fastq.gz"
[3] "trimm_demult/filtN/C_CAGATCCATC_L001_R1_001.fastq.gz"
[4] "trimm_demult/filtN/D_TCGCAATTAC_L001_R1_001.fastq.gz"
[5] "trimm_demult/filtN/E_TTCGAGACGC_L001_R1_001.fastq.gz"
[6] "trimm_demult/filtN/F_TGCCACGAAC_L001_R1_001.fastq.gz"
[7] "trimm_demult/filtN/G_AACCTCATTC_L001_R1_001.fastq.gz"
[8] "trimm_demult/filtN/H_CCTGAGATAC_L001_R1_001.fastq.gz"
[9] "trimm_demult/filtN/I_TTACAACCTC_L001_R1_001.fastq.gz"
[10] "trimm_demult/filtN/J_AACCATCCGC_L001_R1_001.fastq.gz"
[11] "trimm_demult/filtN/K_ATCCGGAATC_L001_R1_001.fastq.gz"
[12] "trimm_demult/filtN/L_TCGACCACTC_L001_R1_001.fastq.gz"
[13] "trimm_demult/filtN/M_CGAGGTTATC_L001_R1_001.fastq.gz"
[14] "trimm_demult/filtN/N_TCCAAGCTGC_L001_R1_001.fastq.gz"
[15] "trimm_demult/filtN/O_TCTTACACAC_L001_R1_001.fastq.gz"
[16] "trimm_demult/filtN/P_TTCTCATTGAAC_L001_R1_001.fastq.gz"
[17] "trimm_demult/filtN/Q_TCGCATCGTTC_L001_R1_001.fastq.gz"
[18] "trimm_demult/filtN/R_TAAGCCATTGTC_L001_R1_001.fastq.gz"
Ahora utilizaremos la función filterAndTrim para eliminar las Ns de las secuencias, y especificamos los path a los reads demultiplexados sin barcodes (guardados en fnFs), y el path a los reads sin Ns que se crearán.
filterAndTrim(fnFs, fnFs.filtN, maxN = 0, multithread = TRUE)
Creating output directory: trimm_demult/filtN
maxN = 0 indica que el número máximo de Ns en cada secuencia es 0
multithreads = TRUE indica que se pueden utilizar varios nucleos en el proceso.
Ahora se eliminarán los adaptadores de los reads. Para esto, primero cargamos la secuencia del adaptador a al objeto FWD para luego obtener la secuencia en todas las orientaciones posibles: directa, inversa, complementaria e inversa complementaria, aplicando una función que llamaremos allOrients:
FWD <- "GATCTTGGTCATTTAGAGGAAGTA" # secuencia del adaptador
allOrients <- function(primer) {
# Create all orientations of the input sequence
require(Biostrings)
dna <- DNAString(primer) # The Biostrings works w/ DNAString objects rather than character vectors
orients <- c(Forward = dna, Complement = complement(dna), Reverse = reverse(dna),
RevComp = reverseComplement(dna))
return(sapply(orients, toString)) # Convert back to character vector
}
FWD.orients <- allOrients(FWD)
FWD.orients
Forward Complement Reverse
"GATCTTGGTCATTTAGAGGAAGTA" "CTAGAACCAGTAAATCTCCTTCAT" "ATGAAGGAGATTTACTGGTTCTAG"
RevComp
"TACTTCCTCTAAATGACCAAGATC"
Ahora contaremos cuántas veces se encuentran las secuencias de los adaptadores en los reads, aplicando una función que llamaremos primerHits.
primerHits <- function(primer, fn) {
# Counts number of reads in which the primer is found
nhits <- vcountPattern(primer, sread(readFastq(fn)), fixed = FALSE)
return(sum(nhits > 0))
}
rbind(FWD.ForwardReads = sapply(FWD.orients, primerHits, fn = fnFs.filtN[[1]]))
Forward Complement Reverse RevComp
FWD.ForwardReads 191594 0 0 0
Se obtuvo una tabla que indica cuantas veces se encontraron las secuencias del adaptador en cada orientación posible.
Para efectivamente eliminar las secuencias adaptadoras, utilizaremos cutadapt, que debe estar instalado en el sistema. Como luego deberemos indicar la ubicación del ejecutable de cutadapt, lo guardaremos en un objeto:
cutadapt <- "/usr/local/bin/cutadapt"
Para ver si podemos eecutar cutadapt podemos correr el siguiente comando, que si no da error significa que el ejecutable de cutadapt funciona correctamente:
system2(cutadapt, args = "--version")
Crearemos un directorio donde se guardarán los reads sin adaptadores:
path.cut <- "demulti/cutadapt"
dir.create(path.cut)
cutFs <- sort(list.files(path.cut, pattern = "_1.fastq.gz", full.names = TRUE))
Finalmente, para eliminar los adaptadores de los reads utilizaremos:
R1.flags <- paste("-g", FWD, "-m 10")
for(i in seq_along(fnFs)) {
system2(cutadapt, args = c(R1.flags, "-n", 2, # -n 2 required to remove FWD and REV from reads
"-o", fnFs.cut[i], # output files
fnFs.filtN[i])) # input files
}
-m 10 indica que se descarten reads menores a 10 bases
Para saber si quedaron reads sin eliminar adaptadores hacemos:
rbind(FWD.ForwardReads = sapply(FWD.orients, primerHits, fn = fnFs.cut[[1]]))
Forward Complement Reverse RevComp
FWD.ForwardReads 0 0 0 0
Todas las columnas con 0 indica que no se detectaron adaptadores en los reads.
Antes de realizar el análisis del metagenoma debemos asegurarnos de trabajar con datos de buena calidad. Para esto, utilizaremos la función filterAndTrim, aplicándola sobre los reads filtrados (sin adaptadores) obtenidos en el paso anterior. Primero se generarán los nombres de los archivos trimmeados y luego se aplica la función. Para los datos de IonTorrent es necesario agregar el parámetro trimLeft=15, que elimina las primeras 15 bases de los reads.
Paso: trimming de los reads por calidad.
Requerimientos:
Se obtiene:
reads filtrados por calidad:
trimLeft = 15: se eliminan las primeras 15 basesmaxLen = 400: se eliminan las bases más allá de la posición 400maxN = 0: máximo de N en los reads es 0compress = TRUE: se comprimirán los archivos trimmeadosPara determinar los parámetros que utilizaremos en el filtrado debemos observar los gráficos de calidad de reads por posición.
# Formato del nombre de los reads
cutFs <- sort(list.files(path.cut, pattern = "_001.fastq.gz", full.names = TRUE))
get.sample.name <- function(fname) strsplit(basename(fname), "_")[[1]][1]
sample.names <- unname(sapply(cutFs, get.sample.name))
head(sample.names)
plotQualityProfile(cutFs)
En el plot se observa que en general la calidad de los reads baja a partir de la posición 400. Según la curva roja, la proporción de reads mayores a esta longitud es baja, por lo que no perderíamos mucha información al trimmear estas bases.
filtFs <- file.path(path.cut, "filtered", basename(cutFs))
out <- filterAndTrim(fwd = cutFs, # archivos para filtrar
filt = filtFs, # nombre de los filtrados
trimLeft = 15, # específico para IonTorrent: elimina las primeras 15 bases de cada read
maxN = 0, # No puede haber Ns
minLen = 50, # Longitud mínima de los reads filtrados
maxLen = 400, # Longitud máxima de los reads filtrados
compress = TRUE, multithread = TRUE,
maxEE = 2, verbose = T)
plotQualityProfile(filtFs)
Plot obtenido a partir de los reads posterior al filtrado de calidad. En este caso la calidad no se ha visto demasiado afectada luego del filtrado realizado. Podrían utilizarse parámetros más robustos, como por ejemplo, eliminar las bases más allá de la posición 400.
A continuación se muestra una tabla que muestra el número de reads previo y post filtrado:
| reads.in | reads.out | |
|---|---|---|
| A_AGGCAATTGC_L001_R1_001.fastq.gz | 197989 | 137116 |
| B_TTAGTCGGAC_L001_R1_001.fastq.gz | 164671 | 127883 |
| C_CAGATCCATC_L001_R1_001.fastq.gz | 188388 | 139982 |
| D_TCGCAATTAC_L001_R1_001.fastq.gz | 329482 | 144849 |
| E_TTCGAGACGC_L001_R1_001.fastq.gz | 240860 | 183122 |
| F_TGCCACGAAC_L001_R1_001.fastq.gz | 281451 | 199745 |
| G_AACCTCATTC_L001_R1_001.fastq.gz | 239688 | 161239 |
| H_CCTGAGATAC_L001_R1_001.fastq.gz | 153471 | 123224 |
| I_TTACAACCTC_L001_R1_001.fastq.gz | 376469 | 273191 |
| J_AACCATCCGC_L001_R1_001.fastq.gz | 276075 | 215029 |
| K_ATCCGGAATC_L001_R1_001.fastq.gz | 213569 | 172837 |
| L_TCGACCACTC_L001_R1_001.fastq.gz | 248175 | 167810 |
| M_CGAGGTTATC_L001_R1_001.fastq.gz | 243272 | 194246 |
| N_TCCAAGCTGC_L001_R1_001.fastq.gz | 243203 | 195955 |
| O_TCTTACACAC_L001_R1_001.fastq.gz | 236879 | 145163 |
| P_TTCTCATTGAAC_L001_R1_001.fastq.gz | 300565 | 237256 |
| Q_TCGCATCGTTC_L001_R1_001.fastq.gz | 147184 | 113679 |
| R_TAAGCCATTGTC_L001_R1_001.fastq.gz | 261204 | 206026 |
dada2 tiene una función para aprender los errores de los reads. plotErrors es una función que permite ver los errores esperados en los reads y los observados.
errF <- learnErrors(filtFs2, multithread = TRUE, verbose = 2)
plotErrors(errF, nominalQ = TRUE) # No quedan como en la figura del Tutorial
Los resultados esperados para este punto se muestran a continuación:
plotErrors: resultados esperados (Segun el Tutorial)
Vamos a eliminar todos los reads repetidos, que aportan información redundante. Este proceso se llama “desreplicar”. Para esto ejecutamos:
derepFs <- derepFastq(filtFs, verbose = TRUE)
names(derepFs) <- sample.names
En este paso es donde dada2 hace su magia. Los detalles del algoritmo pueden consultarse en el paper del paquete: https://www.nature.com/articles/nmeth.3869#methods.
La función dada se aplica sobre los reads limpios, filtrados y trimmeados sin duplicados.
Paso: eliminación de ruido, inferencia de variantes
Requiere:
learnErrors()Se obtiene: una lista con tantos elementos como muestras, donde cada elemento es un objeto de clase dada que contiene los resultados de aplicar la función dada
dadaFs <- dada(derepFs,
err = errF,
multithread = TRUE,
HOMOPOLYMER_GAP_PENALTY = -1,
BAND_SIZE = 32)
multithread = TRUE indica que se utulizarán varios núcleos en paralelo.
HOMOPOLYMER_GAP_PENALTY = -1 está recomendado para datos de Ion Torrent. Se refiere al costo de los gaps en regiones homopoliméricas de más de 3 bases repetidas. Si se deja como NULL estos gaps se tratan igual a los normales
BAND_SIZE = 32 está recomendado para datos de Ion Torrent. Cuando se setea se realiza alineamientos pareados globales (Needleman-Wunsch) por bandas. Las bandas restringen el número cumulativo neto de inserciones en una secuencua contra la otra. El valor por defecto es 16, pero cuando se aplica a datos de pirosecuenciación, con altas tasas de incorporación de indels, debe incrementarse.
Esta función construye una tabla de secuencias similar a la OTU table, que contiene una fila para cada muestra y una columna para cada secuencia única para todas las muestras. Los nombres de las columnas contienen las secuencias encontradas.
seqtab <- makeSequenceTable(dadaFs)
dim(seqtab)
En este paso se eliminan las quimeras.
seqtab.nochim <- removeBimeraDenovo(seqtab,
method = "consensus",
multithread = TRUE, verbose = TRUE)
Al elegir el método consensus, las muestras en la tabla de secuencias son chequeadas independientemente, y se construye una desición consensuada para cada variante de las secuencias.
This function implements a table-specific version of de novo bimera detection. In short, bimeric sequences are flagged on a sample-by-sample basis. Then, a vote is performed for each sequence across all samples in which it appeared. If the sequence is flagged in a sufficiently high fraction of samples, it is identified as a bimera. A logical vector is returned, with an entry for each sequence in the table indicating whether it was identified as bimeric by this consensus procedure.
En este paso, que puede realizarse en cualquier punto del pipeline, vemos cómo ha variado el número de reads a traves de los diferentes pasos del pipeline. Se produce una tabla con tantas filas como muestras y con las columnas correspondientes a los pasos del pipeline: filtrado, denoise y remoción de quimeras.
track %>%
as.data.frame() %>%
mutate(porcent_filt = 100-(filtered/input*100)) %>%
pull(porcent_filt) %>%
summary()
Min. 1st Qu. Median Mean 3rd Qu. Max.
19.07 21.08 23.37 26.92 30.32 56.04
Vemos que en el filtrado se descarta la mayor proporción de los reads, entre 19 y 56% de los reads de partida.
Para la asignación taxonómica utilizaremos la base de datos UNITE, obtenida desde https://unite.ut.ee. Esto se logra con la función assignTaxonomy():
Requiere:
Se obtiene:
head(taxa.print)
Kingdom Phylum Class Order
[1,] "k__Fungi" "p__Mortierellomycota" "c__Mortierellomycetes" "o__Mortierellales"
[2,] "k__Fungi" "p__Ascomycota" "c__Dothideomycetes" "o__Pleosporales"
[3,] "k__Fungi" "p__Ascomycota" "c__Sordariomycetes" "o__Hypocreales"
[4,] "k__Fungi" "p__Ascomycota" "c__Sordariomycetes" "o__Coniochaetales"
[5,] "k__Fungi" "p__Basidiomycota" "c__Tremellomycetes" "o__Filobasidiales"
[6,] "k__Fungi" "p__Ascomycota" "c__Sordariomycetes" "o__Hypocreales"
Family Genus Species
[1,] "f__Mortierellaceae" "g__Mortierella" "s__elongata"
[2,] "f__Didymellaceae" "g__Epicoccum" "s__thailandicum"
[3,] "f__Nectriaceae" "g__Fusarium" NA
[4,] "f__Coniochaetaceae" NA NA
[5,] "f__Piskurozymaceae" "g__Solicoccozyma" "s__phenolica"
[6,] "f__Nectriaceae" "g__Fusarium" "s__kyushuense"
head(taxa.print)
Antes de continuar, veremos el proceso que hemos transitado:
Utilizaremos funciones implentadas en el paquete phyloseq para obtener las tablas necesarias para visualizar los resultados del metagenoma utilizando Shaman: https://shaman.pasteur.fr
Alternativamente puede utilizarse phyloseq para continuar con los análisis y obtener visualizaciones del metagenoma, pero Shaman es más sencillo e intuitivo.
Antes de obtener la OTU table, vamos a construir un data frame que contiene los metadatos de cada muestra:
samples.out <- rownames(seqtab.nochim)
samdf <- read.table("mappingfiledaniellaits_corrected-2.txt", header = T)
samdf <- samdf[1:18,c(1,5,6)]
rownames(samdf) <- samples.out
Ahora, estamos en condiciones de crear un objeto phyloseq que contiene la OTU table, los metadatos y la taxa table:
ps <- phyloseq(otu_table(seqtab.nochim, taxa_are_rows=FALSE),
sample_data(samdf),
tax_table(taxa))
ps
Para construir la taxa table y la OTU table utilizamos las funciones tax_table() y otu_table() respectivamente a partir del objeto phyloseq.
otutable <- otu_table(ps)
otu.table_df <- as.data.frame(otu.table)
rownames(otu.table_df) <- sample.names
taxatable <- tax_table(ps)
taxa.table_df <- as.data.frame(taxatable)
Para exportar las tablas para visualar los resultados en Shaman utilizamos la función write.table sobre las tablas producidas en el paso anterior:
write.table(t(otu.table_df), file = "otu_table.tsv", sep = "\t", quote = F)
write.table(taxa.table_df, "taxa_table.tsv", sep = "\t", quote = F)
Shaman también requiere una tabla de metadatos para realizar las comparaciones entre muestras. Para esto debemos eliminar de la tabla de metadatos todos las variables que contengan información redundante, que contengan datos que se obtienen a partir de transformaciones lineales de otras variables, o que no se utilicen para contrastar las muestras. En este caso elimino las columnas LinkerPrimerSequence, BarcodeSequence, Tiempo y Description, y exporto la tabla resultante:
mf <- read.table("mappingfiledaniellaits_corrected-2.txt", header = T) %>%
select(-LinkerPrimerSequence, -BarcodeSequence, -Tiempo, -Description)
write.table(mf, sep = "\t", "mf.tsv")
Aquí finaliza este pipeline, donde se obtuvieron:
Estas tablas, exportadas como archivos de texto con columnas separadas por tabulaciones puede ser importada en Shaman y analizarse y visualizarse los metagenomas de las muestras, además de realizar comparaciones.